El crecimiento poblacional y la industrialización han llevado al hombre a ocupar prácticamente todos los biomas, impactando gravemente en la diversidad de los sistemas naturales. La biodiversidad es esencial para el mantenimiento de las funciones y servicios ecosistémicos, de los cuales nos beneficiamos, pero a la vez son indispensables para el mantenimiento de la homeóstasis terrestre (Cardinale et al. (2012)). Los esfuerzos para la conservación de estas funciones han llevado a analizar las causas de pérdida de la biodiversidad, siendo las especies invasoras la segunda mayor amenaza a la biodiversidad, después de la pérdida de hábitat (Simberloff et al. (2013)). Para luchar contra esta causa necesitamos políticas de gestión efectivas, pero el esfuerzo de erradicación aumenta con el tiempo de invasión, y por este motivo es imprescindible tratar de trabajar en la previsión y primeras fases de invasión (Van Wilgen, Davies, and Richardson (2014)).
Según la Unión Internacional para la Conservación de la naturaleza (UICN), las especies invasoras son:
“Especie exótica invasora”: especie exótica que se establece en un ecosistema o hábitat natural o seminatural; es un agente de cambio y amenaza la diversidad biológica nativa (UICN, 2000).
Existen diversas hipótesis que tratan de explicar los motivos por los cuales algunas especies se convierten en invasoras. La hipótesis mayoritaria es que la falta de enemigos naturales en el nuevo hábitat dónde la invasora llega, le permite crecer sin control invirtiendo menos en compuestos y estrategias defensivas. Esta es la Enemy Release Hypothesis (Catford, Jansson, and Nilsson (2009)).
Otros autores postulan la regla del 10. De diez especies exóticas que llegan a un nuevo terreno, una se convierte en invasora. Esta hipótesis se basa en el azar y el hecho de que algunas comunidades son mas fáciles de ser invadidas que otras (Richardson and Pyšek (2006)).
Una tercera corriente postula que hay ciertos rasgos que proporcionan a las especies invasoras mayores posibilidades para serlo. Por ejemplo, grandes inversiones en dispersión de semillas, ciclos de vida cortos o crecimiento clonal son algunos de los rasgos propuestos para definir un “buen invasor” (Van Kleunen et al. (2010)).
Ninguna de estas hipótesis es excluyente, por lo que los factores determinantes del éxito invasor podría ser una combinación de aspectos de las tres corrientes (Catford, Jansson, and Nilsson (2009)).
Pese a la divergencia de hipótesis que explican la expansión de las especies invasoras, encontramos una mayor unificación con respecto a las teorías que explican el establecimiento de las especies fuera de su rango nativo. No todas las especies pueden ser encontradas en todos los hábitats. Cada especie tiene sus requerimientos y por lo tanto es posible definir su nicho ecológico en base a las condiciones donde una especie puede vivir, crecer y reproducirse (Grinnell (1917)). El nicho de una especie se define como una región (un hipervolumen de n dimensiones) en un espacio multidimensional basado en factores ambientales que afectan al bienestar de la especie (Hutchinson 1957). La similitud climática entre diferentes países podría ser clave para entender el establecimiento de las especies fuera de su rango nativo y una herramienta muy útil para la valoración del riesgo invasor. Un mejor conocimiento de los factores que determinan si una especie podría convertirse en invasora, sería vital para las políticas de conservación de la biodiversidad, pues la prevención es la acción más económica para erradicar una especie invasora (Kumschick and Richardson (2013)). Saber hasta qué punto la similitud con el marco climático de una región, en nuestro caso la Península Íbérica, puede favorecer a la invasión biológica nos puede permitir definir medidas de prevención mayores.
El objetivo principal de este trabajo es utilizar diferentes herramientas de gestión y análisis de datos para valorar el riesgo invasor en la Península Ibérica, en base a la hipótesis de la similitud climática como principal factor determinante para el establecimiento de espécies exóticas. Los objetivos específicos son:
Marco climatico Ibérico: Espacio multidimensional definido por los valores de diferentes variables bioclimáticas georreferenciados dentro de la Península Ibérica.
Variable bioclimática: Variable descriptora ambiental considerada determinante para el crecimiento, desarrollo y reproducción de la biosfera.
Nicho climático: Espacio multidimensional definido por los valores de diferentes variables bioclimáticas de los puntos de presencia de una especie.
La versión de R con la que se ha trabajado es:
Rv <- R.Version()
Rv$version.string
## [1] "R version 3.5.0 (2018-04-23)"
Antes de empezar, instalaremos los paquetes necesarios a lo largo del documento.
list.of.packages <- c("DataExplorer", "raster", "dismo", "plotly", "fmsb", "factoextra", "MASS","hypervolume","alphahull", "dplyr", "rgdal", "leaflet", "RColorBrewer", "devtools")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
Una vez disponibles todos los paquetes, podemos cargar las respectivas librerías.
library(raster)
library(sp)
library(dismo)
library(ggplot2)
library(plotly)
library(MASS)
require(MASS)
library(fmsb)
library(devtools)
install_github("ggbiplot","vqv")
library(ggbiplot)
library(factoextra)
library(readr)
library(hypervolume)
library(alphahull)
library(DataExplorer)
En primer lugar descargaremos el dataset de la base de datos globales WorldClim que ofrece 19 variables bioclimáticas con una resolución geográfica de 10’ (minutos de grado) ya que resoluciones mayores conllevan demasiado esfuerzo computacional. El argumento var = "bio" de la función getData() reclama dichas variables bioclimáticas. Las variables referentes a la temperatura están multiplicadas por 10.
worldata <- getData("worldclim",var="bio", res=10)
worldata_25 <-getData("worldclim",var="bio",res=2.5)
plot (worldata_25$bio1,
main = "Mapa de la temperatura media (BIO1) global",
xlab = "Longitud (ºC)",
ylab = "Latitud (ºC)")
Los datos de worldata se encuentran en formato rasterizado. Ahora debemos seleccionar la extensión correspondiente a la Península Ibérica. Pasaremos la ventana de latitudes y longitudes que enmarcan la Península Ibérica a la función extent(): (-10º W, 4ºE) y (36ºN, 44ºN). Seguidamente, recortaremos la extensión e sobre los datos worldata utilizando la función crop().
e <- extent(-10, 4, 36, 44)
s.crop <- crop(worldata_25, e)
plot (s.crop$bio1,
main = "Mapa de la temperatura media (BIO1) en la Península Ibérica",
xlab = "Longitud (ºC)",
ylab = "Latitud (ºC)")
El archivo
s.crop también se encuentra en formato rasterizado. Para trabajar más cómodamente, construiremos un dataframe llamado marcoclim, definiendo los nombres de las variables bioclimáticas:
s.crop <- crop(worldata, e)
marcoclim <- as.data.frame(s.crop)
BIOs <- paste('BIO', (1:19), sep = '')
names(marcoclim) <- BIOs
Finalmente, exportamos el dataframe marcoclim en formato *.csv.
write.csv(marcoclim, file="Datos/MarcoClim.csv")
Antes de continuar trabajando, necesitamos hacernos una idea de cómo son los datos, para ello realizamos un pequeño análisis exploratorio y data wrangling.
En la mayoría de ocasiones los datos con los que se trabaja necesitan de una limpieza y una preparación previa para ser analizados. En nuestro caso, eliminamos los missing values del marcoclim.
marcoclim <- marcoclim[complete.cases(marcoclim),] # Elimina NA's.
marcoclim <- na.omit(marcoclim)
head(marcoclim)
## BIO1 BIO2 BIO3 BIO4 BIO5 BIO6 BIO7 BIO8 BIO9 BIO10 BIO11 BIO12 BIO13
## 52 134 86 40 4578 247 35 212 85 175 193 76 1287 157
## 53 132 89 40 4698 248 30 218 81 174 192 72 1259 153
## 54 129 94 41 4862 251 23 228 76 173 191 67 1212 147
## 55 129 99 41 5005 255 19 236 73 192 192 64 1145 137
## 56 129 103 42 5164 260 16 244 71 194 194 62 1073 127
## 57 127 107 42 5311 264 11 253 68 195 195 59 1015 118
## BIO14 BIO15 BIO16 BIO17 BIO18 BIO19
## 52 64 23 422 238 274 394
## 53 63 23 412 235 268 385
## 54 61 22 394 229 258 370
## 55 58 21 368 220 220 346
## 56 55 20 341 209 209 320
## 57 54 19 317 202 202 298
summary(marcoclim) #resumen de cada variable
## BIO1 BIO2 BIO3 BIO4
## Min. : 24.0 Min. : 54.0 Min. :26.00 Min. :2680
## 1st Qu.:117.0 1st Qu.: 91.0 1st Qu.:37.00 1st Qu.:5178
## Median :136.0 Median :102.0 Median :38.00 Median :5746
## Mean :134.6 Mean :100.5 Mean :37.79 Mean :5632
## 3rd Qu.:157.0 3rd Qu.:112.0 3rd Qu.:39.00 3rd Qu.:6236
## Max. :185.0 Max. :129.0 Max. :47.00 Max. :7222
## BIO5 BIO6 BIO7 BIO8
## Min. :150.0 Min. :-76.00 Min. :127 Min. : -1.0
## 1st Qu.:258.0 1st Qu.: -1.00 1st Qu.:237 1st Qu.: 83.0
## Median :285.0 Median : 17.00 Median :266 Median :106.0
## Mean :284.2 Mean : 21.17 Mean :263 Mean :105.8
## 3rd Qu.:312.0 3rd Qu.: 46.00 3rd Qu.:293 3rd Qu.:126.0
## Max. :361.0 Max. : 87.00 Max. :334 Max. :203.0
## BIO9 BIO10 BIO11 BIO12
## Min. :-20.0 Min. : 94.0 Min. :-38.00 Min. : 221.0
## 1st Qu.:181.0 1st Qu.:189.0 1st Qu.: 45.00 1st Qu.: 478.0
## Median :205.0 Median :210.0 Median : 64.00 Median : 593.0
## Mean :193.2 Mean :209.6 Mean : 66.02 Mean : 670.0
## 3rd Qu.:232.0 3rd Qu.:234.0 3rd Qu.: 90.00 3rd Qu.: 800.8
## Max. :275.0 Max. :277.0 Max. :129.00 Max. :1522.0
## BIO13 BIO14 BIO15 BIO16
## Min. : 31.00 Min. : 0.00 Min. :11.00 Min. : 85.0
## 1st Qu.: 62.00 1st Qu.: 5.00 1st Qu.:27.00 1st Qu.:164.0
## Median : 80.00 Median :14.00 Median :37.00 Median :217.0
## Mean : 88.25 Mean :19.62 Mean :39.13 Mean :240.3
## 3rd Qu.:106.00 3rd Qu.:30.00 3rd Qu.:52.00 3rd Qu.:289.0
## Max. :236.00 Max. :90.00 Max. :77.00 Max. :621.0
## BIO17 BIO18 BIO19
## Min. : 9.00 Min. : 16.00 Min. : 66.0
## 1st Qu.: 36.00 1st Qu.: 39.00 1st Qu.:125.0
## Median : 67.50 Median : 76.00 Median :191.0
## Mean : 82.89 Mean : 91.43 Mean :210.4
## 3rd Qu.:117.00 3rd Qu.:130.00 3rd Qu.:262.8
## Max. :304.00 Max. :304.00 Max. :621.0
plot_str(marcoclim) #estructura de los datos
plot_histogram(marcoclim)
cor(marcoclim) # Idea de la correlación de unas variables con otras.
## BIO1 BIO2 BIO3 BIO4 BIO5
## BIO1 1.000000000 -0.003218943 0.050613531 -0.058277439 0.7212966
## BIO2 -0.003218943 1.000000000 0.168482451 0.770769526 0.6318544
## BIO3 0.050613531 0.168482451 1.000000000 -0.477637743 -0.1019689
## BIO4 -0.058277439 0.770769526 -0.477637743 1.000000000 0.5984012
## BIO5 0.721296606 0.631854374 -0.101968935 0.598401199 1.0000000
## BIO6 0.883805142 -0.407885261 0.167118931 -0.492863922 0.3452915
## BIO7 -0.017527499 0.918590949 -0.228272638 0.954158045 0.6709731
## BIO8 0.456888738 -0.108221983 -0.127929509 0.001236178 0.2364828
## BIO9 0.652759550 0.273745339 0.130033272 0.116211839 0.6607126
## BIO10 0.919547558 0.297076397 -0.144429394 0.336239959 0.9167325
## BIO11 0.942102524 -0.264616419 0.191978714 -0.384576525 0.4689499
## BIO12 -0.348841653 -0.559844960 0.182554789 -0.609402384 -0.6593703
## BIO13 -0.100078545 -0.603344167 0.198089465 -0.670610900 -0.5041127
## BIO14 -0.758129653 -0.351421682 -0.020090498 -0.262785988 -0.8207783
## BIO15 0.762416555 -0.015640429 0.099351207 -0.121696887 0.5703048
## BIO16 -0.072684836 -0.575272942 0.231543534 -0.668608154 -0.4659488
## BIO17 -0.731581881 -0.375530061 -0.004061334 -0.293794281 -0.8213257
## BIO18 -0.687875511 -0.448629819 -0.051967481 -0.326059433 -0.8322554
## BIO19 -0.001603785 -0.467961092 0.278680322 -0.608561397 -0.3417890
...
plot_correlation(marcoclim) #para visualizar mejor la correlación
Ahora, podemos representar los datos obtenidos, dependiendo de las variables que nos interesa mostrar:
Para generar un gráfico con temperatura y precipitación:
plot(marcoclim$BIO1~marcoclim$BIO12, pch=19, col="blue", xlab="Precipitación Anual (mm)", ylab="Temperatura anual media (ºC x10)",main="Temperatura y Precipitación en la Península Ibérica")
Podemos graficar las variables teniendo en cuenta su densidad estimada en base a métodos como Kernel Density Estimation. La estimación de la función de densidad para la temperatura media anual seria:
d <- density(marcoclim$BIO1) # devuelve los datos de densidad
plot(d)
Para generar un gráfico de densidad con dos variables usamos la función kde2d(kernel density two dimensions) del paquete MASS, usando las variables de temperatura media anual y precimitación media anual:
dens <- kde2d(marcoclim$BIO1, marcoclim$BIO12, n=300, lims = c(0,200,0,1400))
image(dens)
filled.contour(dens,
color.palette=colorRampPalette(c('gray94','blue','yellow','red','darkred')),
xlab= "Annual Mean Temperature (ºC)",
ylab= "Annual Precipitation (mm)", main = "Marco Climático Ibérico")
Al trabajar con 19 variables, algunas producto lineal de otras como hemos visto en el gráfico de correlación, estamos trabajando con información redundante. Podemos basarnos en el criterio VIF (factores de inflación de varianza), un enfoque simple para identificar la colinealidad entre las variables explicativas para descartar algunas y probar el modelo más adelante con todas las variables y con solo las seleccionadas.
La cuantificación de los factores de inflación de la varianza y el descarte paso a paso se pueden realizar con una función descrita en el blog beckmw:
vif_func<-function(in_frame,thresh=10,trace=T,...){
library(fmsb)
if(any(!'data.frame' %in% class(in_frame))) in_frame<-data.frame(in_frame)
#get initial vif value for all comparisons of variables
vif_init<-NULL
var_names <- names(in_frame)
for(val in var_names){
regressors <- var_names[-which(var_names == val)]
form <- paste(regressors, collapse = '+')
form_in <- formula(paste(val, '~', form))
vif_init<-rbind(vif_init, c(val, VIF(lm(form_in, data = in_frame, ...))))
}
vif_max<-max(as.numeric(vif_init[,2]), na.rm = TRUE)
if(vif_max < thresh){
if(trace==T){ #print output of each iteration
prmatrix(vif_init,collab=c('var','vif'),rowlab=rep('',nrow(vif_init)),quote=F)
cat('\n')
cat(paste('All variables have VIF < ', thresh,', max VIF ',round(vif_max,2), sep=''),'\n\n')
}
return(var_names)
}
else{
in_dat<-in_frame
#backwards selection of explanatory variables, stops when all VIF values are below 'thresh'
while(vif_max >= thresh){
vif_vals<-NULL
var_names <- names(in_dat)
for(val in var_names){
regressors <- var_names[-which(var_names == val)]
form <- paste(regressors, collapse = '+')
form_in <- formula(paste(val, '~', form))
vif_add<-VIF(lm(form_in, data = in_dat, ...))
vif_vals<-rbind(vif_vals,c(val,vif_add))
}
max_row<-which(vif_vals[,2] == max(as.numeric(vif_vals[,2]), na.rm = TRUE))[1]
vif_max<-as.numeric(vif_vals[max_row,2])
if(vif_max<thresh) break
if(trace==T){ #print output of each iteration
prmatrix(vif_vals,collab=c('var','vif'),rowlab=rep('',nrow(vif_vals)),quote=F)
cat('\n')
cat('removed: ',vif_vals[max_row,1],vif_max,'\n\n')
flush.console()
}
in_dat<-in_dat[,!names(in_dat) %in% vif_vals[max_row,1]]
}
return(names(in_dat))
}
}
Aplicamos la función a nuestros datos. Ésta va calculando los VIF, y en cada cálculo quita el que tiene mayor valor hasta quedarse con valores de VIF por debajo del umbral definido (en nuestro caso hasta quedarnos con valores inferiores a 5, un umbral recomendado).
vif_func(in_frame=marcoclim,thresh=5,trace=T)
## var vif
## BIO1 813.522988309366
## BIO2 239.44782413085
## BIO3 34.7762393382378
## BIO4 1092.17479361358
## BIO5 Inf
## BIO6 Inf
## BIO7 Inf
## BIO8 3.96369033534115
## BIO9 4.21583510546391
## BIO10 3011.52623004908
## BIO11 1663.61726148902
## BIO12 236.137816470082
## BIO13 99.181388951579
## BIO14 105.495934592052
## BIO15 25.7859393870437
## BIO16 324.243130786138
## BIO17 172.888788198011
## BIO18 119.650658542968
## BIO19 168.424151605773
##
## removed: BIO5 Inf
...
## [1] "BIO2" "BIO3" "BIO8" "BIO9" "BIO13" "BIO14"
Estas son las variables con las que nos quedariamos para reducir la colinealidad.
Para describir de forma reducida el MCI realizamos una reducción de dimensiones mediante análisis de componentes principales.
Realizamos la PCA con nuestros datos utilizando la función prcomp que R tiene incorporada.
pcaclim <- prcomp(marcoclim, center = TRUE, scale. = TRUE)
print(pcaclim)
## Standard deviations (1, .., p=19):
## [1] 2.932734e+00 2.389085e+00 1.445999e+00 1.073491e+00 8.027324e-01
## [6] 6.757993e-01 4.704042e-01 2.338774e-01 1.637139e-01 1.226069e-01
## [11] 1.002573e-01 7.772701e-02 6.887137e-02 5.791342e-02 4.710744e-02
## [16] 4.006634e-02 3.086699e-02 1.353002e-02 1.181938e-15
##
## Rotation (n x k) = (19 x 19):
## PC1 PC2 PC3 PC4 PC5
## BIO1 0.26057151 0.238012832 -0.13443828 0.036857434 0.24466095
## BIO2 0.17820964 -0.248355819 0.33374275 -0.262346960 0.22369365
## BIO3 -0.02885208 0.127133703 0.15301532 -0.849036098 0.09461437
## BIO4 0.16752742 -0.312197451 0.18058120 0.297613338 0.15748431
## BIO5 0.32210187 -0.013221687 0.13112561 0.082911866 0.26497258
## BIO6 0.14467891 0.355591401 -0.18936869 -0.002411807 0.09563225
## BIO7 0.18968588 -0.293422195 0.27336908 0.080156006 0.17451898
## BIO8 0.12949820 0.003907452 -0.55982629 -0.045122752 0.17893624
## BIO9 0.22447097 0.144032253 0.20457806 -0.053542475 0.34923933
## BIO10 0.31312982 0.103748692 -0.05145848 0.156047032 0.27156093
## BIO11 0.18679938 0.324969551 -0.18207586 -0.052433109 0.14884430
## BIO12 -0.27141600 0.202084486 0.19080726 0.116922878 0.21182655
## BIO13 -0.20471851 0.293152290 0.19246920 0.149316359 0.07935544
## BIO14 -0.30956271 -0.109145795 -0.10041250 0.001946151 0.31115272
## BIO15 0.20969829 0.279824140 0.16464927 0.099991960 -0.30461907
## BIO16 -0.19302030 0.303187985 0.22790285 0.129339110 0.05195185
## BIO17 -0.31266003 -0.087197203 -0.08955716 0.008136523 0.36966961
## BIO18 -0.30985703 -0.073151928 -0.16220531 0.048298047 0.33412666
## BIO19 -0.15581732 0.310362018 0.31426410 0.099878118 0.06984306
...
plot(pcaclim, type = "l")
summary(pcaclim)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 2.9327 2.3891 1.4460 1.07349 0.80273 0.67580
## Proportion of Variance 0.4527 0.3004 0.1100 0.06065 0.03391 0.02404
## Cumulative Proportion 0.4527 0.7531 0.8631 0.92379 0.95770 0.98174
## PC7 PC8 PC9 PC10 PC11 PC12
## Standard deviation 0.47040 0.23388 0.16371 0.12261 0.10026 0.07773
## Proportion of Variance 0.01165 0.00288 0.00141 0.00079 0.00053 0.00032
## Cumulative Proportion 0.99339 0.99626 0.99767 0.99847 0.99899 0.99931
## PC13 PC14 PC15 PC16 PC17 PC18
## Standard deviation 0.06887 0.05791 0.04711 0.04007 0.03087 0.01353
## Proportion of Variance 0.00025 0.00018 0.00012 0.00008 0.00005 0.00001
## Cumulative Proportion 0.99956 0.99974 0.99986 0.99994 0.99999 1.00000
## PC19
## Standard deviation 1.182e-15
## Proportion of Variance 0.000e+00
## Cumulative Proportion 1.000e+00
Obtenemos la desviación estándar y la proporción de la varianza de cada una de los componentes principales basándose en las 19 variables climáticas. Con PC1 y PC2 se explica un 75,44% de la variabilidad. Teniendo en cuenta el objetivo del trabajo, consideramos adecuado quedarnos con los dos primeros componentes.
Además, podemos calcular el peso de cada variable bioclimática para cada componente principal, con pca$rotation
pcaclim$rotation[,1] # PC1
## BIO1 BIO2 BIO3 BIO4 BIO5 BIO6
## 0.26057151 0.17820964 -0.02885208 0.16752742 0.32210187 0.14467891
## BIO7 BIO8 BIO9 BIO10 BIO11 BIO12
## 0.18968588 0.12949820 0.22447097 0.31312982 0.18679938 -0.27141600
## BIO13 BIO14 BIO15 BIO16 BIO17 BIO18
## -0.20471851 -0.30956271 0.20969829 -0.19302030 -0.31266003 -0.30985703
## BIO19
## -0.15581732
pcaclim$rotation[,2] # PC2
## BIO1 BIO2 BIO3 BIO4 BIO5
## 0.238012832 -0.248355819 0.127133703 -0.312197451 -0.013221687
## BIO6 BIO7 BIO8 BIO9 BIO10
## 0.355591401 -0.293422195 0.003907452 0.144032253 0.103748692
## BIO11 BIO12 BIO13 BIO14 BIO15
## 0.324969551 0.202084486 0.293152290 -0.109145795 0.279824140
## BIO16 BIO17 BIO18 BIO19
## 0.303187985 -0.087197203 -0.073151928 0.310362018
Para facilitar el entendimiento de la contribución de cada una de las variables en el análisis de componentes principales, Podemos generar plots con las funciones: -fviz_pca_var -ggbiplot
fviz_pca_var(pcaclim,
col.var = "contrib", # Color by contributions to the PC
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE # Avoid text overlapping
)
g <- ggbiplot(pcaclim, obs.scale = 1, var.scale = 1,
ellipse = TRUE, alpha=0.02, circle = FALSE)
print(g)
También podríamos realizar una PCA solo con las variables que VIF < 5(opcional)
marcoclimVIF <- marcoclim[,c(2,3,8,9,13,15)]
pcaclim2 <- prcomp(marcoclimVIF, center = TRUE, scale. = TRUE)
print(pcaclim2)
## Standard deviations (1, .., p=6):
## [1] 1.3290614 1.2863099 1.0979518 0.9198756 0.5886380 0.4252514
##
## Rotation (n x k) = (6 x 6):
## PC1 PC2 PC3 PC4 PC5
## BIO2 0.59210269 -0.1444876 0.45372154 -0.10819177 0.28852275
## BIO3 0.08306651 0.3467811 0.47623371 0.77388516 -0.02840042
## BIO8 0.16892244 -0.3077534 -0.61995095 0.60118470 -0.01768585
## BIO9 0.50365198 0.4416262 -0.17743188 -0.15289635 -0.70233760
## BIO13 -0.53763457 0.4944007 0.01922602 0.01303984 -0.02774002
## BIO15 0.26695192 0.5698143 -0.38877671 -0.06652759 0.64929780
## PC6
## BIO2 0.57247155
## BIO3 -0.21526599
## BIO8 0.36149418
## BIO9 0.05624511
## BIO13 0.68206148
## BIO15 -0.16397395
plot(pcaclim2, type = "l")
summary(pcaclim2)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 1.3291 1.2863 1.0980 0.9199 0.58864 0.42525
## Proportion of Variance 0.2944 0.2758 0.2009 0.1410 0.05775 0.03014
## Cumulative Proportion 0.2944 0.5702 0.7711 0.9121 0.96986 1.00000
Como podemos comprobar, usando el criterio VIF la variabilidad explicada en el PCA con dos componentes resulta del 58.05%, un porcentaje inferior que cuando no eliminamos las variables redundantes, por tanto decidimos no trabajar siguiendo este criterio.
Llegados a este punto, pretendemos evaluar si las condiciones climáticas de los diferentes países del mundo se parecen a las de la Península Ibérica. Así se podrá determinar después qué especies podrían encontrar condiciones favorables en nuestro país actuando como invasoras y qué especies españolas podrían serlo para otros países. Con estas evidencias se controlarían mejor la entrada y salida de especies invasoras, mejorando las políticas de comercio y transporte por ejemplo. Por eso proyectamos los datos del marco climático de estos países sobre el nuestro.
Preparamos los datos de la Península Ibérica de los dos primeros componentes:
marcoclim <- marcoclim[complete.cases(marcoclim),] # Elimina NA's
e <- pcaclim$x[,1:2]
edf <- as.data.frame(e)
edf$pais <- "PenIberica" #Añadimos una varaible de País, con todas las etiquetes "PenIberica"
head(edf)
## PC1 PC2 pais
## 52 -5.632903 2.4546673 PenIberica
## 53 -5.482798 2.0945529 PenIberica
## 54 -5.206638 1.5848376 PenIberica
## 55 -4.514175 1.1294438 PenIberica
## 56 -3.987905 0.6762998 PenIberica
## 57 -3.618972 0.1304936 PenIberica
En este blog, encontramos como se realiza la descarga y recorte de los países. Para llamar a un país concreto se usa la nomenclatura de tres letras definida aquí. Hay una base de datos con los nombres y codigos de los 241 países que incluye la lista. Podemos importar esos datos:
CountryCodes <- read_delim("CountryCodes.csv",
";", escape_double = FALSE, trim_ws = TRUE)
## Parsed with column specification:
## cols(
## Country = col_character(),
## A2 = col_character(),
## A3 = col_character(),
## `Number,` = col_number()
## )
#A3 son los codigos para descargar el país.
Descargamos los límites políticos de un país, por ejemplo Dinamarca, de la siguiente forma:
denmark <- getData('GADM', country='DNK', level=0)
s<- crop(worldata, denmark)
s2 <- mask(s, denmark)
s2<-as.matrix(s2)
marcoclim.denmark <- as.data.frame(s2)
names(marcoclim.denmark) <- BIOs
marcoclim.denmark <- marcoclim.denmark[complete.cases(marcoclim.denmark),]
h<- crop(worldata_25, denmark)
plot (h$bio1,
main = "Mapa de la temperatura media (BIO1) en Dinamarca",
xlab = "Longitud (ºC)",
ylab = "Latitud (ºC)")
Así, obtenemos un
dataframe con las variables bioclimáticas de este país.
Proyectamos estos valores sobre el PCA del marco climático ibérico. Generamos la variable país, con el nombre de la etiqueta del mismo. Además combinamos los dos data frames, con los datos de PC1 y PC2 de España y el país seleccionado.
i <- predict (pcaclim, newdata = marcoclim.denmark)
i <- i[,1:2]
idf <- as.data.frame(i)
idf$pais <- "Dinamarca"
df <- rbind(edf,idf)
Generamos un gráfico para visualizar el grado de solapamiento de ambas nubes de puntos usando la función ggplot
p<-ggplot(df, aes(x=PC1, y=PC2, col=pais)) + geom_point() + ggtitle("Proyección de Dinamarca sobre el PCA del Marco Climático Ibérico")
p + labs(x = "PC1 (46.3%)", y = "PC2 (29.1%)" )
Para entender mejor el solapamiento podemos generar gráficos de densidad:
Para el Marco Climático Ibérico:
dense <- kde2d(edf$PC1, edf$PC2, n=300, lims = c(-12,7,-10,10))
image(dense)
filled.contour(dense,
color.palette=colorRampPalette(c('gray94','blue','yellow','red','darkred')),
xlab= "PC1",
ylab= "PC2", main = "Marco climático Ibérico")
Para Dinamarca:
densi <- kde2d(idf$PC1, idf$PC2, n=300, lims = c(-12,7,-10,10))
image(densi)
filled.contour(densi,
color.palette=colorRampPalette(c('gray94','blue','yellow','red','darkred')),
xlab= "PC1",
ylab= "PC2", main = "Marco climático Ibérico del país proyectado")
Pero, ¿Cómo podemos cuantificar el grado de solapamiento, teniendo en cuenta la densidad de los datos?
Con el fin de determinar el grado de solapamiento climático de un país cualquiera respecto al MCI, creamos hipervolúmenes a partir de los dos primeros componentes principales, que explican un 75.44% de la varianza, y los comparamos, obteniendo así el volumen compartido y otros índices de similitud.
Para ello, utilizaremos el paquete hypervolume (Benjamin, Christine, et al. (2017))(Benjamin, Babich, et al. (2017)), un paquete diseñado para estimar la forma y volumen de bases de datos de más de una dimension y realizar operaciones con ellos: intersección, unión, determinación de los componentes únicos, tests de inclusión y detección de agujeros. Usa una aproximación basada en geometria estocástica ([CRAN] (https://cran.r-project.org/web/packages/hypervolume/hypervolume.pdf)).
Con el objetivo de determinar valores de solapamiento del MCI con todos los países del mundo, crearemos un bucle para contrastar país a país. En primer lugar definiremos el hipervolumen definido por los puntos en el PCA del MCI, a lo que denominaremos marco estático. El marco dinámico estará constituido por aquellos elementos en bucle, dinámicos, de los cuales iremos sacando un conjunto de outputs. En las diferentes pestañas definiremos estos elementos.
Una vez obtenidos los datos del MCI, realizado el análisis de componentes principales y creado un data frame con los resultados de éste, definiremos el primer hipervolumen. El hipervolumen de n-dimensiones es una medida que generaliza el concepto de volumen a espacios de dimensión. En nuestro caso utilizaremos un Gaussian KDE(kernel density estimation) en el que todos los puntos contribuyen a la densidad de la probabilidad global ya que testeando vemos que es el mejor algoritmo de inclusion de puntos por nuestra tipologia de distribuciones.
En la función hypervolume_gaussian, por defecto, se le asigna un umbral del 95%, asegurando la máxima probabilidad en la estimación de la densidad. Además generaremos un plot para visualizar el hipervolumen en el que se genera un centroide (en rojo al centro del volumen).
hv1sp = hypervolume_gaussian(data=subset(df,pais=="PenIberica")[1:2730, 1:2])
##
## Building tree...
## done.
## Ball query...
##
## done.
##
## Building tree...
## done.
## Ball query...
##
## done.
## Requested probability quantile 0.950000, obtained 0.949316 - setting threshold value 0.000001.
## For a closer match, you can increase num.thresholds in hypervolume_threshold.
plot(hv1sp,show.contour=TRUE,contour.kde.level=0.01)
El siguiente paso será generar un hipervolumen para cada uno de los países que testearemos frente al MCI. Por ejemplo, creamos un hipervolúmen con los datos de Dinamarca.
hv2dn = hypervolume_gaussian(data=subset(df,pais=="Dinamarca")[1:238, 1:2])
##
## Building tree...
## done.
## Ball query...
##
## done.
##
## Building tree...
## done.
## Ball query...
##
## done.
##
## Building tree...
## done.
## Ball query...
##
## done.
##
## Building tree...
## done.
## Ball query...
##
## done.
## Requested probability quantile 0.950000, obtained 0.949491 - setting threshold value 0.000183.
## For a closer match, you can increase num.thresholds in hypervolume_threshold.
plot(hv2dn,show.contour=TRUE,contour.kde.level=0.01)
Una vez obtenidos los dos hipervolúmenes, utilizaremos la función hypervolume_set y hypervolume_overlap_statistics. La primera genera valores estadísticos computando la intersección, la unión y los componentes únicos(diferencia) de los hipervolúmenes, que podemos observar llamando la función get_volume
La segunda función calcula los métricos de superposición para los dos hipervolúmenes. Imprime directamente los índices de similitud entre los dos volúmenes:
Similitud de Jaccard: Volumen de la intersección de un país respeto al otro, dividido por el volumen de la unión de estos países.
Similitud de Sorensen: Dos veces el volumen de la interesección de un país sobre el otro, dividido por la suma del volumen de ambos países.
Fracción Única 1: Volumen del componente único del primer país dividido por el volumen total de éste.
Fracción única 2: Volumen del componente único del segundo país dividido por el volumen total de éste.
Además calcularemos la distancia entre centroides usando la función hypervolume_distance
setH <-hypervolume_set(hv1sp, hv2dn, check.memory = FALSE)
## Choosing num.points.max=25955 (use a larger value for more accuracy.)
## Using minimum density of 159.367907
## Retaining 15092 points in hv1 and 796 points in hv2.
## Beginning ball queries...
##
## Building tree...
## done.
## Ball query...
##
## done.
##
## Building tree...
## done.
## Ball query...
##
## done.
## Finished ball queries.
volumes<-get_volume(setH)
overl<- hypervolume_overlap_statistics(setH)
dist<- hypervolume_distance(hv1sp, hv2dn, type='centroid')
Con los datos estadísticos obtenidos anteriormente, podemos generar una lista para facilitar la visualización de los diferentes índices y volúmenes.
values <- c(volumes, overl, dist)
En esta figura ilustramos los diferentes índices obtenidos:
Métricos de contraste de similitud
En este punto, implementaremos un sencillo for loop para testear la similitud del marco climático de cada país respecto al de la Península Ibérica de forma rápida y sistemática. El MCI será un valor fijo en este bucle, mientras que los inputs seran los datos climáticos de los diferentes países a testear. A partir de este bucle, generaremos como output un DataFrame que recogerá los volúmenes e índices respecto al MCI de la mayoría de países del mundo. Primero cargaremos una lista con los nombres y los acrónimos de todos los países del mundo:
CountryCodes <- read.csv("CountryCodes.csv", header=FALSE, sep=";", skip = 1)
Decimos “mayoría” de países porque los datos climáticos de algunos de ellos no presentan una distribución geográfica suficientemente amplia como para ser procesados correctamente por la función hypervolume_gaussian del paquete hypervolume. Para evitar que el bucle se interrumpa, debemos eliminar manualmente dichos países:
black_list = list('ASM', 'AIA', 'ATA', 'ABW', 'BMU', 'BVT', 'IOT', 'CYM', 'CXR', 'CCK', 'GIB', 'GRD', 'KIR', 'LIE', 'MAC', 'MDV', 'MLT', 'MHL', 'FSM', 'MCO', 'MSR', 'NRU', 'ANT', 'NIU', 'NFK', 'MNP', 'PCN', 'RUS', 'SHN', 'KNA', 'LCA', 'SPM', 'VCT', 'SMR', 'SCG', 'SYC', 'TKL', 'TUV', 'UMI', 'VAT', 'VGB', 'WLF')
CountryCodesSum <- CountryCodes[ ! CountryCodes$V3 %in% black_list, ]
Como puede verse, la mayoría de los países de la lista black_list son muy pequeños, como por ejemplo SMR (San Marino), GIB (Gibraltar) o VAT (Ciudad del Vaticano). También aparecen muchos estados insulares, como es el caso de CYM (Islas Caimán), BMU (Bermudas) o SYC (Seychelles). ATA (Antártida) aparece en esta lista, no por tener una superficie pequeña (que no es el caso), sino por la escasa distribución geogràfica de sus datos climáticos.
Un caso aparte es el de RUS (Rusia), que ha sido colocada en la black_list debido a la extremadamente amplia distrubución geográfica de sus datos climáticos (debido a su gran extensión). Los ordenadores con los que hemos trabajado no disponían de memoria suficiente como para cargar los datos generados por la función hypervolume_gaussian aplicada a Rusia.
A continuación mostraremos el código del blucle que hemos implementado. Es importante aclarar que el siguiente bloque de código ha sido configurado para no ser ejecutado mediante la opción eval = FALSE:
El bucle recorre todos los países en un lapso de unas 7 horas aproximadamente. Esto trabajando con una workstation de la siguientes caracteríticas:
Como ya hemos avanzado, el output de nuestro bucle es un DataFrame llamado Volumes que incluye todos los índices y volúmenes proporcionados por la función hypervolume_gaussian. Finalmete, editaremos brevemente el DataFrame Volumes con la intención de llegar a un DataFrame un poco más leíble y friendly.
CountryCodes$V2 <- NULL
CountryCodes$V4 <- NULL
Volumes <- read.csv("Volumes.csv", header = T)
Volumes$X <- NULL
Vol.Names <- c('Vol.MCI', 'Vol.Country', 'Vol.Intersection', 'Vol.Union', 'Vol.Unique.MCI', 'Vol.Unique.Country', 'Jaccard', 'Sorensen', 'NormVol.MCI', 'NormVol.Country', 'Dist.Centroid', 'Country.Acronym')
names(Volumes) <- Vol.Names
rownames(Volumes) <- Volumes$Country.Acronym
names(CountryCodes) = c('Country.Name', 'V3')
Volume_dataset <- merge(Volumes, CountryCodes, by.x = "Country.Acronym", by.y = "V3")
write.csv(Volume_dataset, file = "Volumes_dataset.csv")
Los diferentes índices y valores obtenidos anteriormente nos permiten ordenar los países de forma que podamos visualizar aquellos que presentan un mayor grado de similitud climática.
Cada índice nos revela características distintas de la intersección entre los volúmenes. Vamos a explorar los resultados obtenidos considerando el top 10 de países por cada índice.
Volumes_dataset <- read.csv("Volumes_dataset.csv")
rVolCount <- arrange(Volumes_dataset, desc(Vol.Country)) # Ordenado de mayor volumen a menor
rVolCount <- rVolCount$Country.Acronym # Guardamos la columna con los nombres de los países
rintersection <- arrange(Volumes_dataset, desc(Vol.Intersection)) #Ordenado de mayor interesección a menor
rintersection <- rintersection$Country.Acronym
rvoluniqmci <- arrange(Volumes_dataset, Vol.Unique.MCI) # Ordenado de menor a mayor volumen único en la Península Iberica.
rvoluniqmci <- rvoluniqmci$Country.Acronym
rvoluniqcount <- arrange(Volumes_dataset, Vol.Unique.Country) # Ordenado de menor a mayor volumen único en la Península Ibérica.
rvoluniqcount <- rvoluniqcount$Country.Acronym
rjaccard <- arrange(Volumes_dataset, desc(Jaccard)) # Ordenado de mayor a menor grado de similitud de Jacacard
rjaccard <- rjaccard$Country.Acronym
rsorensen <- arrange(Volumes_dataset, desc(Sorensen)) # Ordenado de mayor a menor grado de similitud de Sorensen
rsorensen <- rsorensen$Country.Acronym
rnormvolmci <- arrange(Volumes_dataset, NormVol.MCI) # Ordenado de menor a mayor fracción única de MCI
rnormvolmci <- rnormvolmci$Country.Acronym
rnormvolcount <- arrange(Volumes_dataset, NormVol.Country) # Ordenado de menor a mayor fracción única del país contrastado
rnormvolcount <- rnormvolcount$Country.Acronym
rdistance <- arrange(Volumes_dataset, Dist.Centroid) # Ordenado de menor a mayor distancia entre centroides
rdistance <- rdistance$Country.Acronym
ranking <- data.frame(rVolCount, rintersection, rvoluniqmci, rvoluniqcount, rjaccard, rsorensen, rnormvolmci, rnormvolcount, rdistance)
names(ranking) <- c('Vol.Country', 'Vol.Intersection', 'Vol.Unique.MCI', 'Vol.Unique.Country', 'Jaccard', 'Sorensen', 'NormVol.MCI', 'NormVol.Country', 'Dist.Centroid')
topten <- ranking[1:10,] #Sacamos los 10 primeros países.
print(topten)
## Vol.Country Vol.Intersection Vol.Unique.MCI Vol.Unique.Country Jaccard
## 1 IND CHL CHL TCA ESP
## 2 MMR ESP ESP NLD ARG
## 3 COL USA USA MYT GRC
## 4 BTN ARG ARG URY ZAF
## 5 CHL BOL BOL ALA ITA
## 6 PNG AUS AUS PLW LBN
## 7 ECU MEX MEX BHR PRT
## 8 USA GRC GRC AND ALB
## 9 CRI ZAF ZAF LUX TUR
## 10 TWN LBN LBN BRB AUS
## Sorensen NormVol.MCI NormVol.Country Dist.Centroid
## 1 ESP CHL URY ESP
## 2 ARG ESP NLD ARG
## 3 GRC USA ESP GRC
## 4 ZAF ARG PRT LSO
## 5 ITA BOL LSO PRT
## 6 LBN AUS GRC ZAF
## 7 PRT MEX BEL ALB
## 8 ALB GRC ALB URY
## 9 TUR ZAF LUX ITA
## 10 AUS LBN FRA LBN
La mejor forma de presentar los valores obtenidos de similitud entre países es diseñar un mapa mundial con diferentes colores para aquellas zonas con mayor similitud climática con la Península Ibérica, y por lo tanto mayor riesgo invasor.
Para ello, usaremos el paquete leaflet, que nos permite la creación de mapas interactivos. En primer lugar descargamos los datos de fronteras y polígonos de todos los países del mundo desde el portal Web de Thematicmapping.
Necesitaremos también el paquete rgdal, para obtener csv con los datos descargados, con la función readOGR. Los datos deben estar en un repositorio real del ordenador para cargarlos, por ejemplo así:
#world_spdf=readOGR(dsn= "getwd() ",layer="TM_WORLD_BORDERS_SIMPL-0.3")
world_spdf=readOGR(dsn= "C:/WorldBorders",layer="TM_WORLD_BORDERS_SIMPL-0.3")
Combinamos los datos de los polígonos con la variable que queremos visualizar:
ÍNDICE DE JACCARD:
world_spdf@data$number<-c(1:246)
dfJac <- Volumes[,c(7,12)]
names(dfJac) <- c("Jaccard", "ISO3")
world_spdf@data<-merge(world_spdf@data, dfJac, by= "ISO3", all.x= TRUE)
world_spdf@data<-arrange(world_spdf@data, number)
world_spdf@data$Jaccard[ which(world_spdf@data$Jaccard == 0)] = NA
pal <- colorQuantile("YlOrRd", NULL, n = 9)
state_popup <- paste0("<strong> País: </strong>",
world_spdf$NAME,
"<br><strong>Indice de Jaccard: </strong>",
world_spdf$Jaccard)
leaflet(data = world_spdf) %>%
addProviderTiles("CartoDB.Positron") %>%
addPolygons(fillColor = ~pal(Jaccard),
fillOpacity = 0.8,
color = "#BDBDC3",
weight = 1,
popup = state_popup) %>%
addLegend(pal= pal, values= ~Jaccard, opacity=0.9, title = "Jaccard (%)", position = "bottomleft" )
El resultado de éste codigo genera un gráfico interactivo que permite visualizar con zoom los valores del índice en los diferentes países. Aquí mostramos una imagen PNG del mismo, pero para acceder al gráfico interactivo seguir este enlace. Los valores NA corresponden a una similitud nula.
Grado de similitud climatica respecto la Península IbéricaSimiltud climática en base al índice de Jaccard
La Ley 42/2007, de 13 de diciembre, del Patrimonio Natural y de la Biodiversidad creó, en su artículo 61.1, el Catálogo Español de Especies Exóticas Invasoras en el que se incluyen todas aquellas especies exóticas invasoras que constituyan, de hecho, o puedan llegar a constituir una amenaza grave para las especies autóctonas, los hábitats o los ecosistemas, la agronomía, o para los recursos económicos asociados al uso del patrimonio natural (Ministerio de Agricultura , Pesca y Alimentación).
El Catálogo Español de Especies Exóticas Invasoras contempla 181 especies de diferentes grupos taxonómicos provinientes de diferentes países.
Con el objetivo de comprobar si aquellos países con mayor coincidencia climática son fuente de más especies actualmente invasoras, se ha realizado una búsqueda en los portales de ISSG y GBIF para definir los países de origen de las actuales invasoras.
El resultado de esta búsqueda es el siguiente data frame:
origeninvasoras <- read_excel("origeninvasoras.xlsx")
Recontamos el numero de ocurrencias de cada país y lo mezclamos con los datos de similitud:
countriesorigen <-origeninvasoras[,7:56]
b<-as.matrix(countriesorigen)
a<-as.vector(b)
a.t<-table(a)
a.df <- as.data.frame(a.t)
df <- mutate_all(a.df, funs(toupper))
names(df) <- c("Country.Name", "NumSpecies")
withsp<-merge(df, Volume_dataset, by="Country.Name", all.y = TRUE)
withsp$NumSpecies <- as.integer(withsp$NumSpecies)
Seleccionamos el top 10 de países que han aportado mayor número de especies invasoras:
ttorigin <- arrange(withsp, desc(NumSpecies))
ttorigin <- data.frame(ttorigin$Country.Acronym, ttorigin$NumSpecies)
ttorigin[1:10,] #Sacamos los 10 primeros países.
## ttorigin.Country.Acronym ttorigin.NumSpecies
## 1 USA 45
## 2 CHN 29
## 3 MEX 26
## 4 BRA 24
## 5 ARG 22
## 6 JPN 21
## 7 CAN 19
## 8 PRY 19
## 9 URY 18
## 10 IND 17
Podemos representar en forma de correlación los resultados obtenidos, para determinar si la mayor similitud climática explica el número de especies invasoras que aporta ese país:
Distáncia entre centroides vs. Número de especies invasoras aportadas
withsp <- withsp[complete.cases(withsp),]
p <- ggplot(withsp, aes(x=NumSpecies, y=Dist.Centroid)) +
geom_point() +
geom_smooth(method =lm, se=FALSE )
p + labs(x = "Número de especies invasoras en España del país considerado", y = "Distáncia entre centroides climáticos" )
Valor de Índice de Jaccard vs. Número de especies invasoras aportadas
p <-ggplot(withsp, aes(x=NumSpecies, y=Jaccard)) +
geom_point() +
geom_smooth(method =lm, se=FALSE)
p + labs(x = "Número de especies invasoras en España del país considerado", y = "Índice de Jaccard" )
A parte de la similitud del marco climático, hay otros factores que pueden resultar determinantes para la determinación del riesgo invasor. La cantidad de mercancías y personas que se desplazan de un país a otro, podrían ser clave para entender en que medida realmente podrían ocurrir los desplazamientos de especies exóticas. Además, en los resultados presentados en el apartado anterior, hay una variabilidad muy grande, probablemente por trabajar con datos de especies invasoras de diferentes grupos taxonómicos. Todo esto requiere la obtención de nuevos datos y desarrollo de nuevos modelos. De forma exploratoria presentaremos algunos modelos simples explorando datos disponibles como son la biodiversidad de cada país (que implicaría mayor probabilidad de existir una espécie que pueda ser invasora) y la distáncia entre el país y la Península Ibérica (lo que probablemente facilite también la llegada de especies alóctonas).
BIODIVERSIDAD: Existen datos de biodiversidad de especies normalizados por el tamaño del país: el National Biodiversity Index. Los valores estan normalizados entre 1 (Màximo: Indonesia) y 0 (Mínimo: Groenlándia).
Importamos los datos
CountryData <- read_excel("Datos/CountryDataSum.xlsx")
names(CountryData)[3] <- "Country.Acronym"
complete<- merge(Volume_dataset, CountryData, by = "Country.Acronym")
complete <-merge(withsp, complete, by="Country.Acronym")
Generamos el gráfico:
complete$nbi <- as.numeric(complete$NBI)
## Warning: NAs introducidos por coerción
p <- ggplot(complete, aes(x=NumSpecies, y=nbi)) +
geom_point() +
geom_smooth(method = loess, se= TRUE)
p + labs(x = "Número de especies invasoras en España del país considerado", y = "Índice de Biodiversidad nacional" )
## Warning: Removed 25 rows containing non-finite values (stat_smooth).
## Warning: Removed 25 rows containing missing values (geom_point).
Los 10 países con mayor valor de biodiversidad son:
nbirank <- arrange(complete, desc(nbi))
nbirank <- data.frame(nbirank$Country.Acronym, nbirank$nbi)
nbirank[1:10,] #Sacamos los 10 primeros países.
## nbirank.Country.Acronym nbirank.nbi
## 1 IDN 1.000
## 2 COL 0.935
## 3 MEX 0.928
## 4 BRA 0.877
## 5 ECU 0.873
## 6 AUS 0.853
## 7 VEN 0.850
## 8 PER 0.843
## 9 CHN 0.839
## 10 CRI 0.820
DISTÁNCIA GEOGRÁFICA: A través del portal Web Geodatos hemos podido obtener los valores de distancia geográfica de todos los países respecto la Península Ibérica.
complete$distance.km <- as.numeric(complete$distance.km)
p <- ggplot(complete, aes(x=NumSpecies, y=distance.km)) +
geom_point() +
geom_smooth(method = lm, se= FALSE)
p + labs(x = "Número de especies invasoras en España del país considerado", y = "Distáncia respecto el país de origen (km)" )
Durante el desarrollo de este trabajo, hemos aprendido a utilizar distintas herramientas computacionales para valorar el grado de similitud climática de la Península Ibérica con el resto del mundo, con el objetivo de estimar el riesgo invasor. Distintos paquetes de Ry algunas bases de datos disponibles on-line, referentes a variables bioclimáticas de las distintas regiones del planeta, nos han permitido desarrollar un estudio que consideramos muy útil para determinar riesgo invasor, y que sienta las bases de posibles proyectos futuros de evaluación de este tipo. Estos resultados, ayudarían al entendimiento de la problemática relacionada con la pérdida de hábitat y diversidad y podrían mejorar políticas de gestión que protejan o controlen la entrada y salida de especies con potencial invasor, protegiendo así la biodiversidad de nuestro país y también del resto.
Fijándonos en los rankings que hemos obtenido, podemos encontrar diferentes países con más riesgo respecto a la Península Ibérica, dependiendo del métrico o índice al que atendamos. En nuestro caso, hemos determinado el índice de Jaccard como el parámetro más representativo de la similitud de volúmenes para la selección de estos países con condiciones más similares y por tanto de riesgo. Como podemos apreciar en el último mapa, lo países ordenados de mayor a menor riesgo para la Península Ibérica son:
El resultado de este ranking tiene sentido pues todos los países listados, contienen clima mediterráneo, o bien comparten cercanía con España, o latitud geográfica.
Teniendo en cuenta los otros rankings de número de especies fuente, de biodiversidad y de distancia geográfica, podemos ir más allá en la evaluación del riesgo invasor. En primer lugar, la mitad de países que aparecen en la lista de los 10 países que más especies invasoras han aportado a España, también aparecen en los top 10 de los parámetros de similitud medidos. Por ejemplo, Argentina es el quinto país del que provienen más especies invasoras y a la vez es el país con mayor índice de Jaccard, es decir, más parecido climáticamente a la península ibérica. En segundo lugar, atendiendo la biodiversidad presente respecto a la superficie del país, podemos observar como uno de los países que se encuentra en el top 10 del índice de Jaccard, Australia, también lo está en el top 10 de mayor biodiversidad. Es decir, Australia presenta un riesgo añadido a parte de la similitud climática ya que aumenta la probabilidad de que alguna de las especies llegue a ser invasora en la península ibérica. Otro aspecto interesante a comentar en referencia a los datos de biodiversidad, sería que de la lista del top 10 de países con mayor volumen climático (proyectado sobre la península ibérica), algunos se corresponden con los países de mayor biodiversidad mundial, confirmando la relación ya descrita de mayor biodiversidad en lugares climáticamente más diversos.
Si asumimos que el número de especies invasoras en España viene definido por distintos factores, podríamos evaluar la correlación de cada uno de ellos con ese número. Por ejemplo, según nuestros resultados observamos que a mayor similitud (índice de Jaccard y distancia entre centroides), mayor número de especies invasoras provenientes de ese país. A mayor biodiversidad nacional se ha observado una tendencia a incrementar el número de especies invasoras que han llegado a la península ibérica. Por último, atendiendo a la distancia entre otros países con España observamos que a menor distancia, menor número de especies invasoras seguramente porque los países o regiones cercanas forman parte de la distribución ecológica de las especies que se encuentran en nuestro país o dicho de otra forma, no tendrían carácter invasor por su probable cercanía filogenética y similitud biogeográfica.
Finalmente, nos gustaría resaltar las múltiples posibilidades en cuanto a las variables que se podrían tener en cuenta para seguir evaluando el riesgo, teniendo en cuenta lo desarrollado y obtenido en este trabajo. Por ejemplo, testeo del riesgo invasor de especies animales o plantas a nivel individual y taxonómico o la determinación del riesgo en base al movimiento de ciudadanos o mercancía entre países.
Visto en perspectiva, el presente trabajo representa una primera aproximación de un método para valorar el riesgo invasor desde un punto de vista del data science, una ciencia emergente con gran cantidad de herramientas y en constante evolución que con el tiempo, permitirán el desarrollo de modelos más transversales y de mayor complejidad.
Benjamin, Blonder, Morrow Cecina Babich, Maitner Brian, Harris David J., Lamanna Christine, Violle Cyrille, Enquist Brian J., and Kerkhoff Andrew J. 2017. “New Approaches for Delineating N-Dimensional Hypervolumes.” Methods in Ecology and Evolution 9 (2): 305–19. doi:10.1111/2041-210X.12865.
Benjamin, Blonder, Lamanna Christine, Violle Cyrille, and Enquist Brian J. 2017. “Using N-Dimensional Hypervolumes for Species Distribution Modelling: A Response to Qiao et Al. ().” Global Ecology and Biogeography 26 (9): 1071–5. doi:10.1111/geb.12611.
Cardinale, B.J., J.E. Duffy, A. Gonzalez, D.U. Hooper, C. Perrings, P. Venail, A. Narwani, et al. 2012. “Biodiversity Loss and Its Impact on Humanity.” Nature 486 (7401): 59–67. doi:10.1038/nature11148.
Catford, J.A., R. Jansson, and C. Nilsson. 2009. “Reducing Redundancy in Invasion Ecology by Integrating Hypotheses into a Single Theoretical Framework.” Diversity and Distributions 15 (1): 22–40. doi:10.1111/j.1472-4642.2008.00521.x.
Grinnell, Joseph. 1917. “The Niche-Relationships of the California Thrasher.” The Auk 34 (4). American Ornithological Society: 427–33. http://www.jstor.org/stable/4072271.
Kumschick, S., and D.M. Richardson. 2013. “Species-Based Risk Assessments for Biological Invasions: Advances and Challenges.” Diversity and Distributions 19 (9): 1095–1105. doi:10.1111/ddi.12110.
Richardson, D.M., and P. Pyšek. 2006. “Plant Invasions: Merging the Concepts of Species Invasiveness and Community Invasibility.” Progress in Physical Geography 30 (3): 409–31. doi:10.1191/0309133306pp490pr.
Simberloff, Daniel, Jean-Louis Martin, Piero Genovesi, Virginie Maris, David A. Wardle, James Aronson, Franck Courchamp, et al. 2013. “Impacts of Biological Invasions: What’s What and the Way Forward.” Trends in Ecology & Evolution 28 (1): 58–66. doi:https://doi.org/10.1016/j.tree.2012.07.013.
Van Kleunen, M., W. Dawson, D. Schlaepfer, J.M. Jeschke, and M. Fischer. 2010. “Are Invaders Different? A Conceptual Framework of Comparative Approaches for Assessing Determinants of Invasiveness.” Ecology Letters 13 (8): 947–58. doi:10.1111/j.1461-0248.2010.01503.x.
Van Wilgen, B.W., S.J. Davies, and D.M. Richardson. 2014. “Invasion Science for Society: A Decade of Contributions from the Centre for Invasion Biology.” South African Journal of Science 110 (7-8). doi:10.1590/sajs.2014/a0074.